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Liquid  water  transport  in  a  polymer  electrolyte  fuel  cell  (PEFC)  is  a  major  issue  for  automotive  applications. 
Mist  flow  with  tiny  droplets  suspended  in  gas  has  been  commonly  assumed  for  channel  flow  while  two- 
phase  flow  has  been  modeled  in  other  cell  components.  However,  experimental  studies  have  found  that 
two-phase  flow  in  the  channels  has  a  profound  effect  on  PEFC  performance,  stability  and  durability. 
Therefore,  a  complete  two-phase  flow  model  is  developed  in  this  work  for  PEFC  including  two-phase 
flow  in  both  anode  and  cathode  channels.  The  model  is  validated  against  experimental  data  of  the  wetted 
area  ratio  and  pressure  drop  in  the  cathode  side.  Due  to  the  intrusion  of  soft  gas  diffusion  layer  (GDL) 
material  in  the  channels,  flow  resistance  is  higher  in  some  channels  than  in  others.  The  resulting  flow 
maldistribution  among  PEFC  channels  is  of  great  concern  because  non-uniform  distributions  of  fuel  and 
oxidizer  result  in  non-uniform  reaction  rates  and  thus  adversely  affect  PEFC  performance  and  durability. 
The  two-phase  flow  maldistribution  among  the  parallel  channels  in  an  operating  PEFC  is  explored  in 
detail. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  portability,  compactness,  zero  emission  and  high  power  out¬ 
put  at  low  temperature  has  made  the  polymer  electrolyte  fuel  cell 
(PEFC)  one  of  the  most  potent  replacements  for  the  internal  com¬ 
bustion  engines  [1].  This  new  focus  has  led  to  an  urgent  need  for 
identification,  understanding,  prediction,  control,  and  optimization 
of  various  transport  and  electrochemical  processes  that  occur  on 
disparate  length  scales  in  fuel  cells  [2].  Recent  studies  have  shown 
that  among  all  transport  phenomena,  the  two-phase  transport  of 
water  in  the  PEFC  to  maintain  water  balance  is  the  most  critical  to 
cell  performance. 

A  typical  PEFC  and  its  components  are  schematically  displayed 
in  Fig.  1  a.  A  PEFC  model  should  consider  transport  phenomena  with 
electrochemical  kinetics  and  charge  transport  of  both  electrons  and 
protons  in  disparate  length  and  time  scales  [2  ].  The  need  for  detailed 
model  validation  has  been  increasingly  acknowledged  by  the  PEFC 
research  community  because  the  global  current-voltage  curve  is 
largely  inadequate  to  differentiate  various  transport  and  electro¬ 
chemical  processes. 

Much  effort  has  been  directed  toward  PEFC  modeling  [2-5]. 
Although  water  is  essential  for  membrane  proton  conductivity, 
excess  water  can  initiate  channel  flooding,  blocking  the  pores  of 
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GDL  and  catalyst  layer  and  hampering  the  reactant  transport  [6,7]. 
Channel  flooding  refers  to  a  situation  where  a  substantial  frac¬ 
tion  of  liquid  water  accumulates  in  gas  channels.  Given  the  low 
startup  temperatures  (room  temperature)  for  automobile  appli¬ 
cations,  two-phase  flow  is  unavoidable  for  automobile  fuel  cells. 
Therefore,  its  understanding  and  prediction  is  critical  for  good 
PEFC  design.  Two-phase  transport  in  a  PEFC  consists  of  three  sub¬ 
processes:  (1)  liquid  water  accumulation  and  transport  in  catalyst 
layers,  (2)  two-phase  transport  in  GDL  and  MPLs,  along  with  the 
interfacial  coverage  at  the  GDL  surface,  and  (3)  two-phase  flow 
in  gas  channels.  All  two-phase  flow  modeling  efforts  in  PEFC  in 
the  literature  were  on  the  first  two  sub-processes  [8-13].  Ample 
experimental  evidence,  however,  indicates  that  the  channel  flood¬ 
ing  plays  a  pivotal  role  in  water  management,  particularly  at  low 
current  densities  where  gas  velocity  is  insufficiently  low  to  drain 
liquid  water  out  of  the  channels.  The  low  load  regime  is  particu¬ 
larly  important  for  PEFC  engines  due  to  its  potential  of  high  energy 
conversion  efficiency  and  it  is  most  frequently  used. 

The  flooding  problem  is  compounded  by  GDL  intrusion  into  the 
channels.  PEFC  is  operated  under  high  compression  to  minimize  the 
contact  resistance  between  the  land  and  GDL.  The  high  compression 
pressure  pushes  the  softer  GDL  material  into  the  channel,  blocking 
the  channels  partially  as  shown  in  Fig.  lb.  Compression  pressure 
is  highest  at  the  edge  (near  the  location  of  tightening  bolts)  and 
therefore  GDL  intrusion  is  most  severe  at  the  edge  channels.  Flow 
through  the  intruded  channel  is  reduced  under  the  same  pressure 
drop,  thereby  making  it  more  difficult  to  flush  liquid  water  out  of 
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Nomenclature 

A  area 

Q  local  concentration  of  species  i,  mol  m-3 

Dj  diffusion  coefficient  of  species  i,  m2  s-1 

d  diameter 

/  current  density,  A  cm-2 

I<  permeability  of  the  porous  media 

k  relative  permeability  of  the  phase 

M  molecular  weight 

mfk  mass  fraction  of  species  k  in  liquid  phase 
n  normal  direction 

P  pressure,  Pa 

5  source  term  in  the  governing  equations 

s  liquid  water  saturation/volume  fraction 

T  absolute  temperature,  I< 

u  velocity 

x  x  coordinate,  |xm 

y  y  coordinate,  f±m 

z  z  coordinate,  p,m 

Greek  letters 

a  net  water  transport  coefficient 

p  density 

s  porosity 

y  correction  factor 

k  electrolyte  conductivity,  S  m-1 

A  relative  mobility 

a  surface  tension  coefficient  (N  nrr1 ) 

§  Stoichiometry  at  gas  channel  inlet  (anode  or  cath¬ 

ode) 

v  kinematic  viscosity 

Subscripts  and  superscripts 
g  gas  phase 

1  liquid  species 

eff  effective 

sat  saturation 

r  relative 

c  convective  correction 

mem  membrane 

in  inlet 

ph  phase 

avg  average 

H20  water 

02  oxygen 


the  channel.  This  exacerbates  channel  flooding  and  accelerates  the 
mass  transport  loss  as  well  as  leads  to  operational  instability.  The 
worst  scenario  is  total  blockage  of  a  channel  by  liquid  water.  The 
maldistribution  of  flow  in  parallel  channels  has  profound  perfor¬ 
mance  and  durability  implications,  and  a  serious  loss  of  efficiency 
is  possible  as  the  whole  channel  is  lost  due  to  blockage.  Hence,  for  a 
PEFC  to  maintain  stable  performance,  flooding  of  the  channels  must 
be  avoided. 

Channel  flooding  in  PEFCs  has  received  increased  attention  in 
the  fuel  cell  community  [14-26].  A  friction  factor  similar  to  that  in 
the  case  of  laminar  flow  through  a  circular  channel  was  proposed 
for  the  PEFC  channels  [14].  Preferential  entry  of  inlet  flow  in  paral¬ 
lel  channels  was  attributed  partly  to  the  formation  of  recirculation 
vortex  at  the  inlet  [15].  Single-phase  flow  in  the  channel  has  been 
used  widely  in  the  gas  channel  for  design  purposes  [16,17,18]  but 
experimental  observations  indicate  presence  of  liquid  water  in  the 
PEFC  gas  channels  [19,20].  Attempts  have  been  made  to  use  ana¬ 


lytically  calculated  dry  length  of  a  gas  channel  as  a  design  tool  for 
PEFC  channel  flow  field  design  [21].  Others  used  the  void-in-fluid 
(VOF)  method  [22-26]  to  compute  two-phase  flow  in  the  gas  chan¬ 
nels.  Unfortunately,  these  models  are  computationally  expensive 
and  have  not  been  integrated  with  the  other  components  of  a  PEFC. 

Few  experimental  investigations  on  flow  distribution  in  the  PEFC 
channels  exist  in  literature  [  14,19,20,27].  While  mist  flow  model  and 
film  flow  model  have  been  used  for  the  extreme  cases  of  high  gas 
velocity  (liquid  volume  fraction  <0.1%)  and  low  gas  velocity/highly 
hydrophilic  channel  wall  (liquid  volume  fraction  >10%),  respec¬ 
tively,  a  general  model  covering  a  common  range  of  liquid  fraction 
and  capable  of  capturing  flow  maldistribution  was  absent.  Such  a 
general  model  will  enable  the  prediction  of  channel  flooding,  two- 
phase  flow  maldistribution  in  multiple,  parallel  channels,  and  the 
flow-field  effect  on  liquid  water  removal  in  operating  PEFCs. 

In  the  present  work,  we  first  couple  a  recently  developed  two- 
phase  channel  flow  model  with  other  two-phase  models  for  the 
catalyst  layer  and  GDL  previously  developed  in  our  laboratory  to 
form  a  complete  two-phase  model  for  an  entire  PEFC.  The  channel 
two-phase  model  [28,29]  is  based  on  the  framework  of  multi¬ 
phase  mixture  model  (M2  model)  and  capable  of  predicting  the 
liquid  water  volume  fraction  and  pressure  along  the  flow  direction. 
Then,  the  complete  model  is  validated  against  experimental  data 
of  wetted  area  ratio  and  pressure  drop  over  a  range  of  operating 
conditions.  Finally,  this  complete  two-phase  model  is  employed  to 
study  the  effects  of  GDL  intrusion  and  manifold  design  on  reducing 
flow  maldistribution. 

2.  Mathematical  model 

A  PEFC  consists  of  seven  sub-regions  -  anode  gas  channel,  anode 
GDL,  anode  catalyst  layer,  ionomeric  membrane,  cathode  catalyst 
layer,  cathode  GDL  and  cathode  gas  channel.  In  addition,  the  elec¬ 
tron  transport  through  bipolar  plates  may  be  important  in  some 
cases  [30].  The  membrane  is  a  solid-state  electrolyte  with  water 
and  proton  co-transport  taking  place  through  its  ionomer  phase. 
Full  descriptions  of  electrochemical  and  transport  phenomena  in  a 
PEFC  already  exist  in  the  literature  [2]  and  are  not  repeated  here. 

Two-phase  flow  and  transport  in  a  PEFC  are  governed  by  the 
laws  of  momentum,  mass,  energy,  species  and  charge  conserva¬ 
tion.  Under  non-isothermal,  two-phase  conditions  the  conservation 
equations  of  mass,  momentum,  energy,  species  and  charge  equa¬ 
tions  in  the  PEFC  can  be  written  as  [30-32] 

Mass:  V  •  (pu)  =  0  (1) 


Momentum:  -^V  •  (puu)  =  -VP  -  V  •  (pr)  +  SU  (2) 


Energy : 

V  •  (yjpCpuT)  =  V  ■  (keffVT)  +  ST 

(3) 

Species  : 

V  ■  (ycuCk)  =  V  •  (Dg  effVCg) 

-V- 

\(<  cg\-l 
\M«  pgJJ\ 

+  S/,. 

(4) 

Charge  (electrons) :  V  •  (creffV<2>s)  +  S^s  =  0 

(5) 

Charge  (protons) :  V  •  (/<effV<2>e)  +  =  0 

(6) 

The  source  terms  are  tabulated  in  Table  1.  Details  about  these 
equations  and  the  source  terms  are  available  in  the  literature 
[31,32  ].  The  present  modeling  approach  is  to  view  all  components  in 
a  PEFC  as  porous  media.  Specially,  we  model  flow  channels  of  typical 
dimension  between  0.2  and  1  mm  as  a  structured  porous  media  or 
a  bundle  of  straight  capillary  tubes.  Hence,  we  apply  the  two-phase 
flow  theory  based  on  extended  Darcy’s  law  to  describe  two-phase 
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(a) 


Cathode  BipolarPlate 


Cathode  Catalyst  Layer 
(02+4H*+4e  ->2H20) 


(b) 


GDL  intruction  in  the  channel 


Land 


Land 


Fig.  1.  (a)  Schematic  of  a  multi-channel  counter-flow  PEFC.  (b)  Schematic  of  a  typical  GDL  intrusion  in  a  channel,  (c)  Sections  of  the  computational  mesh. 


flow  and  transport  throughout  an  entire  PEFC  including  the  flow 
channels.  Furthermore,  for  convenience  of  numerical  implementa¬ 
tion,  the  two-phase  theory  in  porous  media  is  formulated  into  M2 
model  without  making  any  additional  assumptions  [8,11,12,32,33]. 
The  validity  of  this  porous  medium  approach  for  channel  two-phase 
flow  encountered  in  PEFCs  has  been  explored  in  detail  elsewhere 
[28,29]. 

Assuming  a  local  thermodynamic  equilibrium  between  the 
liquid  and  vapor  phases  in  the  two-phase  region,  the  water  concen¬ 


tration  in  the  vapor  phase  could  be  taken  as  equal  to  the  saturation 
water  concentration  that  solely  depends  on  the  temperature.  Since 
the  liquid  phase  consists  of  only  water,  the  water  mole  fraction  in 
the  liquid  is  equal  to  unity.  Therefore,  the  total  water  concentration 
CH 2°  to  be  solved  by  Eq.  (4)  can  be  written  as 


Ch2°  =  sCh2o  +  (1  _  s)ch2o  =  s  ( +  (1  -  s)CsHa^°  (7) 

\Mh2o  J 


Table  1 

Source  terms  for  the  conservation  equations  in  each  sub-region. 


/  Mk  =  chemical  formula  of  species  k 

Note:  Electrochemical  reaction  where  =  ne  (  Sfe  =  stoichiometry  coefficient  .  In  PEFC,  there  are  (anode)  H2-2H+  =  2e_;  (cathode)  2H20-02-4H+  =4e_. 

\n  =  number  of  electrons  transferred 
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Table  2 

Cell  geometry  and  properties. 


Description 

Value 

Cell  length 

100.0  mm 

Gas  channel  depth 

0.5  mm 

Gas  channel  width 

1.0  mm 

Land  width 

1.0  mm 

Anode/cathode  GDL  thickness 

0.20  mm 

Anode/cathode  catalyst  layers  thickness 

0.010  mm 

Porosity  of  anode/cathode  GDL,  sGDL 

0.6 

Porosity  of  anode/cathode  catalyst  layers,  8Gl 

0.6 

Volume  fraction  of  ionomer  in  anode/cathode  catalyst 

0.18 

layers,  se 

Hydraulic  permeability  of  anode/cathode  GDL,  KGDL 

3.0  x  10“14  m2 

Hydraulic  permeability  of  membrane,  /Wm 

5.0  x  10-20  m2 

Contact  angle  of  anode/cathode  channel,  GDL  and  catalyst 

60°,  110°,  110° 

layers,  6 

Contact  resistance  between  catalyst  layer  and  GDL,  RGDL 

1.0  x  10-6  £2  m2 

Anode/cathode  inlet  pressure,  Pin 

2.0  atm 

Cell  temperature,  Tcen 

80°C 

where  s  is  the  liquid  water  saturation  or  volume  fraction.  However, 
this  summation  does  not  indicate  a  mixing  in  molecular  level.  For 
any  two-phase  flow  in  any  control  volume  we  can  calculate  total 
water  concentration  as  in  Eq.  ( 7)  but  it  may  not  be  possible  to  solve  a 
conservation  equation  for  it.  Since  pressure  gradient  and  flow  veloc¬ 
ity  are  linearly  coupled  in  porous  media,  it  is  possible  to  solve  for 
the  total  water  concentration  instead  of  individual  phases  in  porous 
media.  From  Eq.  (7),  the  liquid  saturation  or  volume  fraction  can  be 
calculated  from  [31  ] 


CH2°-CsHa2° 

A/M h2o  -  CsHa2t° 


(8) 


Within  the  M2  model  framework,  the  kinematic  viscosity  of  the 

two-phase  mixture  is  defined  as 


v  = 


fk ri  +  krg\ 

V  yi  y§ ) 


(9) 


where  v\  and  vg  are  the  kinematic  viscosities  of  liquid  and  gas 
phases,  respectively,  while  kr\  and  krg  are  the  relative  permeabilities 
of  liquid  and  vapor  phases,  respectively. 


Table  3 

Simulation  parameters. 


Description 


Value 


Exchange  current  density  x  ratio  of  reaction  surface  to 
catalyst  layer  volume  in  anode  side,  az^ 

Exchange  current  density  x  ratio  of  reaction  surface  to 
catalyst  layer  volume  in  cathode  side,  azjj^ 

Activation  energy  for  oxygen  reduction  reaction  in  cathode 
side,  £a 

Reference  hydrogen  molar  concentration,  cH2,ref 
Reference  oxygen  molar  concentration,  c02!ref 
Anodic  and  cathodic  transfer  coefficients  for  hydrogen 
oxidation  reaction  (HOR) 

Cathodic  transfer  coefficient  for  oxygen  reduction  reaction 
(ORR) 

Dry  membrane  density,  pmem 

Equivalent  weight  of  electrolyte  in  membrane,  EW 

Faraday  constant,  F 

Universal  gas  constant,  Ru 

Surface  tension,  cr 

Liquid  water  density,  p\  (80  °C) 

Liquid  water  viscosity,  /q 

Effective  electronic  conductivity  in  catalyst  layers,  <rG l 
Effective  electronic  conductivity  in  GDL,  crGDL 
Electronic  conductivity  in  current  collector,  crland 
Catalyst  coverage  coefficient,  nc 
Diffusivity  correction  factor,  n 


1.0  x  10* * * 9  Am-3 

2.0  x  104  Am-3 

73, 269 J  mol”1 

40.88  mol  irr3 
40.88  mol  irr3 
aa=ac  =  l 

ac  =  l 

2,000  kg  ur3 
1.1  kg  mol-1 
96,487  C  mol”1 
8.314 J  (molK)-1 
0.0625  Nm-1 
972  kg  m-3 
3.5xlO-4Nsirr2 
1,000  Sm-1 
10,000  Sm-1 
20,000  Sm-1 
2.0 
2.3 


Stoichiometry 


Fig.  2.  Wetted  area  ratio  versus  air  stoichiometry  for  Jav  g  =  0.2  A  cm-2. 


Water  is  the  only  species  that  could  be  present  in  both  liquid 
and  vapor  states.  Therefore,  the  correction  factor  takes  into  account 
contributions  due  to  both  phases  for  water  conservation  equation. 
Other  species  can  be  present  in  only  the  gas  phase.  The  species 
correction  factor  (yc)  can  be  derived  as  follows  [11]: 


P  ( Mg  c  \ 

Ch2°  \M^o  +  PgLsatJ 

Mg 

,  Pg(l  -5) 


(10) 


The  mobility  of  each  phase  A^s)  and  Ag(s)  are  defined  in  Eqs.  (11 ) 
and  (12)  in  terms  of  the  relative  permeabilities  of  the  liquid  and  gas 
and  phases: 


Ms) 


K\  (S)/Vy 

krl{s)/vl  +  krg{s)/vg 


(11) 


Ag(s)  =  1  -  Ms) 


(12) 


Two  most  important  parameters  to  describe  two-phase  flows 
through  a  porous  medium  are  relative  permeability  and  capillary 
pressure  functions.  The  relative  permeabilities  used  in  this  work 


Stoichimetry 


Fig.  3.  Wetted  area  ratio  versus  air  stoichiometry  for  Jav  g  =  0.5  A  cm-2. 


Wetted  area  ratio 
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Stoichiometry 


Fig.  4.  Wetted  area  ratio  versus  air  stoichiometry  for  JaVg  =  0.8  A  cm-2. 


Stoichiometry 


Fig.  5.  Pressure  drop  validation  with  33%  GDL  intrusion. 


(a)  Anode  Outlet/ 
Cathode  Inlet 


(b)  Anode  Outlet/ 
Cathode  Inlet 


(c) 


Anode  Outlet/ 
Cathode  Inlet 
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Fig.  6.  Current  density  (Am  2 )  contours  in  the  membrane. 
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(a)  Cathode  Inlet 
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Cathode  Outlet 


(b)  Cathode  Inlet 
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(c)  Cathode  Inlet 


Fig.  7.  Pressure  (Pa)  contours  at  the  mid-height  of  cathode  channels. 


can  be  expressed  as 

kr[  =  snk  (13) 

krg  =  (l-s)n*  (14) 


where  exponent  nk  varies  depending  on  the  flow  conditions  and 
porous  medium  microstructures.  Different  values  have  been  used 
for  nk  [11,31,34]  depending  on  flow  situations.  More  recently  it  was 
found  [28]  that  in  the  PEFC  flow  channels  the  value  of  nk  equal 
to  5  results  in  the  best  match  with  experimental  pressure  drop 
data.  However,  for  different  channel  geometry  and  dimension  fresh 
calibration  with  experimental  data  may  be  required.  Therefore,  in 
the  present  case,  we  use  nk  equal  to  5.0  in  the  channel  and  4.0 
otherwise. 

Capillary  pressure  is  usually  expressed  as  Leverett  function  of 
the  liquid  saturation  such  that: 

Pc  =  acos(0c)(-J  J(s)  (15) 

We  use  the  Leverett  J(s)  function  given  by  [5] 

...  (  1 .417(1  -  s)  -  2.120(1  -  s)2  +  1 .263(1  -  s)3  for(9c  >  90° 

J(S>  ~  \  1.417s  -  2.120s2  +  1.263s3  for0c  <  90° 

(16) 


The  absolute  permeability  I<  for  the  flow  channels  can  be  com¬ 
puted  by  numerical  experiments  of  simulating  single-phase  flow 
through  the  flow  channels.  In  the  case  of  GDL  intrusion,  the  chan¬ 
nel  cross-sectional  area  decreases.  This  would  result  in  reduction  in 
the  absolute  permeability.  Consider  that  the  absolute  permeability 
through  a  minichannel  can  be  expressed  by  [28] 


I<  = 


(17) 


where  c  is  the  shape  factor  describing  various  cross-section  geome¬ 
tries  of  the  channels  and  dh  is  the  hydraulic  diameter  of  the  channel. 
Thus,  one  can  use  the  following  approximate  expression  for  the 
absolute  permeability  in  the  intruded  channels: 


^intruded  —  ^unintruded(l  ^) 


(18) 


where  8  is  the  fraction  reduction  in  the  channel  height  due  to  GDL 
intrusion.  The  value  of  c  could  be  calculated  analytically  for  only 
a  few  regular  cross-sections,  for  all  other  cross-sections  it  is  cal¬ 
culated  by  matching  the  experimentally  obtained  pressure  drop. 
Therefore,  the  effect  of  GDL  roughness  on  pressure  drop  along  chan¬ 
nel  could  also  be  incorporated  in  this  model. 
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Fig.  8.  Liquid  water  volume  fraction  contours  at  the  mid-height  of  cathode  channels. 


3.  Boundary  conditions 


The  inlet  velocity  (uin)  is  calculated  in  terms  of  cathode  stoi¬ 
chiometry  (§c),  average  current  density  (Javg),  active  area  of  the 
membrane  (Amem),  molar  concentration  of  oxygen  (C°2)  and  the 
cross-sectional  area  of  channels  (Ain)  as  the  following: 


_  ^cfavgAnem 
“in  =  4FC°2/\in 


(19) 


where  F  is  the  Faraday  constant. 

The  molar  concentrations  of  species  at  the  inlet  are  determined 
by  the  inlet  pressure  and  humidity  according  to  the  ideal  gas  law. 
The  exit  boundary  is  assumed  to  be  fully  developed  namely: 


P  =  V  ref, 


In  this  work  the  reference  pressure  is  2  atm. 
At  all  wall  boundaries  we  have: 


(20) 


u.h  =  0, 


(21) 


4.  Numerical  procedures 


The  governing  equations,  Eqs.  (3)-(8),  along  with  their  appro¬ 
priate  boundary  conditions  are  discretized  by  the  finite  volume 
method  and  solved  in  a  commercial  flow  solver,  Fluent  (version 


6.1.22),  by  SIMPLE  algorithm.  The  source  terms  and  physical  prop¬ 
erties  are  implemented  using  the  user-defined  functions  (UDF) 
available  with  commercial  CFD  software  Fluent  6.1.  Overall  species 
balance  and  charge  balance  are  checked  in  addition  to  the  equa¬ 
tion  residuals  as  important  convergence  criteria.  The  cell  geometry 
and  the  simulation  parameters  are  listed  in  Tables  2  and  3.  In 
the  simulations  to  be  presented  below,  all  species  imbalances 
are  less  than  1%  and  residuals  smaller  than  10-5.  Extensive  grid 
independence  test  have  been  performed  by  Meng  and  Wang 
[35]. 

5.  Results  and  discussion 

5  A.  Experimental  validation 

The  most  important  predictive  capabilities  of  our  model  include 
the  fraction  of  the  wet  GDL-channel  interface  and  the  total  two- 
phase  pressure  drop  through  the  flow  channels.  Prediction  of  these 
quantities  is  therefore  validated  against  the  recent  experiments  of 
Hussaini  and  Wang  [27].  A  seven  parallel  channel  14  cm2  active  area 
cell,  shown  schematically  in  Fig.  1 ,  is  used  in  these  experiments.  The 
channels  are  1  mm  x  0.5  mm  in  cross-section  and  100  mm  in  length. 
The  membrane  active  area  is  14  mm  in  width  and  100  mm  in  length 
and  the  membrane  is  30  [xm  thick  Gore  composite  membrane.  The 
computational  mesh  for  this  cell  consists  of  0.15  million  grid  points. 
Sections  of  this  computational  mesh  are  shown  in  Fig.  lc. 
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Fig.  9.  Liquid  water  saturation  contours  at  the  midsection  of  the  cathode  GDL. 


The  wet  area  on  the  GDL-channel  interface  (covered  partially  by 
liquid  water)  is  quantified  from  the  flow  visualization  experiment 
on  the  seven-channel  PEFC.  The  ratio  of  the  wet  area  to  the  total 
GDL-channel  interface  area  is  defined  as  the  “wetted  area  ratio”. 
Here  liquid  saturation  of  0.01  or  greater  is  taken  as  wet  area.  An 
average  of  the  wetted  area  ratio  over  the  entire  cell  including  all 
the  channels  is  taken  as  the  wetted  area  ratio  for  the  PEFC  at  that 
operating  condition.  The  wetted  area  ratio  is  a  strong  function  of 
three  operational  parameters  -  the  average  current  density,  the  air 
stoichiometry  and  the  relative  humidity  at  the  cathode  inlet.  At  low 
inlet  relative  humidity,  ability  of  the  cathode  flow  to  carry  water  in 
vapor  form  is  high.  Therefore,  low  wetted  area  ratio  is  expected 
at  low  relative  humidity.  Similarly,  at  high  air  stoichiometry  flow 
can  carry  more  water  in  vapor  form  and  as  a  result  wetted  area 
ratio  is  expected  to  decrease.  Similarly  at  low  current  density  as 
the  flowrate  is  lower  we  expect  high  wetted  area  ratio.  High  wet¬ 
ted  area  ratio  is  associated  with  high  liquid  volume  fraction  in  the 
cathode  gas  channels  and  therefore  local  surface  coverage  by  liquid 
water  is  expected  to  be  greater.  Absence  of  adequate  experimen¬ 
tal  data  in  this  case  restricts  further  validation  of  our  model.  The 
computed  wetted  area  ratio  in  the  cathode  side  is  compared  with 
the  experimental  data  over  a  range  of  these  operating  conditions 
in  Figs.  2-4.  Fig.  2  depicts  the  model-experimental  comparison  for 
relative  humidities  of  67%,  45%  and  25%  at  /aVg  =  0.2  A  cm-2  with 
cell  temperature  of  80  °C.  Overall,  good  agreement  is  found  (RMS 
(Root  Mean  Square)  error  is  less  than  0.06).  At  a  medium  current 
density  of  0.5  A  cm-2,  the  match  shown  in  Fig.  3  is  also  reason¬ 


able  (RMS  error  is  less  than  0.08)  except  that  the  present  model 
overpredicts  at  relative  humidity  of  67%  for  large  air  stoichiometry. 
The  match  between  experimental  results  and  numerical  predic¬ 
tions  shown  in  Fig.  4  is  again  good  (RMS  error  is  less  than  0.04) 
at  a  high  current  density  of  0.8  A  cm-2.  These  comparisons  shown 
in  Figs.  2-4  demonstrate  that  our  model  can  predict  the  wetted  area 
ratio  reasonably  well  for  a  range  of  practical  operating  conditions. 

The  wetted  area  ratio  less  than  unity,  as  displayed  in  Figs.  2-4, 
physically  implies  that  the  flow  channels  are  partially  dry  (free 
of  any  liquid  water)  and  partially  wet  due  to  liquid  water  accu¬ 
mulation,  and  that  there  exists  a  dry-to-wet  transition  within  the 
channel  length  [32,36].  The  results  shown  in  Figs.  2-4  clearly  indi¬ 
cate  that  the  flow  channels  have  a  longer  wet  portion  under  high 
inlet  humidity,  lower  air  stoichiometry  and  lower  current  density, 
as  expected. 

The  pressure  drop  along  the  cathode  channels  of  the  PEFC 
is  compared  between  experimental  measurements  [27]  and  the 
present  calculations,  as  shown  in  Fig.  5.  Note  that  GDL  intrusion 
to  reduce  cross-sectional  areas  of  flow  channels  must  be  accounted 
for,  as  intrusion  of  soft  GDL  material  into  flow  channels  is  inevitable 
during  assembly  of  PEFCs.  While  it  is  difficult  to  quantify  the  degree 
of  GDL  intrusion  in  situ,  we  found  33%  GDL  intrusion  at  the  edge 
channels  yields  a  reasonable  match  with  the  experimental  pressure 
drop  data  for  all  current  densities,  as  shown  in  Fig.  5.  The  agree¬ 
ment  is  poor,  particularly  for  the  low  current  density  of  0.2  A  cm-2. 
This  may  be  due  to  the  inadequacy  of  the  present  M2  model  to 
describe  significant  two-phase  flow  occurring  at  low  current  den- 
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(a)  Anode  Outlet  (b)  Anode  Outlet 


Fig.  10.  Pressure  (Pa)  contour  at  the  section  of  anode  channels. 


sities.  Future  investigation  is  needed  to  improve  model  accuracy  in 
this  regime. 

5.2.  Flow  maldistribution  and  its  effect  on  cell  performance 

While  experimental  validation  has  been  carried  out  for  system 
parameters,  in  this  subsection  we  reveal  the  two-phase  flow  mald¬ 
istribution  in  an  operating  PEFC  and  the  impact  on  cell  performance 
using  the  present  model.  All  cases  presented  below  use  the  cell 
temperature  of  80  °C  and  relative  humidity  of  67%,  most  typical  of 
automotive  applications.  The  effect  of  GDL  intrusion  (Toray  Paper 
TGPFI-060)  is  studied  for  high  current  density  (0.8  A  cm-2 )  and  low 
stoichiometry  (2.0).  GDL  intrusion  of  17%  and  33%  at  the  two  edge 
channels  are  considered  to  investigate  the  variation  of  flow  maldis¬ 
tribution  and  the  pressure  penalty.  The  area  maldistribution  due  to 
GDL  intrusion  changes  the  distribution  of  liquid  water  among  chan¬ 
nels  and  furthermore,  interacts  with  the  non-linear  characteristics 
of  two-phase  flow  to  result  in  significant  non-uniform  distribution 
of  reactants.  The  results  of  perfectly  symmetric  and  intruded  flow 
channels  are  compared  and  analyzed. 

The  current  density  distributions  at  the  center  of  the  mem¬ 
brane  for  perfect  channels,  17%  intruded  channels  and  33%  intruded 
channels  are  displayed  in  Fig.  6.  The  contour  shows  somewhat  sym¬ 
metric  current  distribution  over  all  seven  channels  for  the  perfect 


channel  case.  The  slight  non-uniformity  among  channels  is  caused 
by  the  fact  that  both  inlet  and  outlet  manifolds  are  included  in 
the  present  computations.  Flow  resistance  in  the  outlet  manifolds 
of  anode  and  cathode  differs  from  that  of  the  inlet  manifold  due 
to  the  presence  of  two-phase  flow  and  changes  with  operating 
conditions.  For  the  17%  intruded  channels,  there  is  an  additional 
reduction  in  the  flow  cross-section  of  the  two  edge  channels  and 
the  local  current  density  there  is  lower.  This  uneven  distribution 
in  current  density  from  channel  to  channel  is  more  severe  for  the 
33%  intrusion.  Clearly,  this  is  caused  by  low  reactant  flows  through 
the  intruded  channels.  Furthermore,  it  is  seen  from  Fig.  6  that  local 
current  density  is  usually  smaller  over  the  area  of  a  land  separating 
two  flow  channels.  Note  that  the  two  edge  channels  are  bounded  by 
a  half  land  on  the  edge  of  the  cell.  The  low  current  density  regions 
resulting  from  flow  maldistribution  decreases  the  cell  performance 
as  well  as  the  utilization  of  precious-metal  catalysts. 

The  pressure  contours  along  the  cathode  channels  are  shown 
in  Fig.  7  for  the  same  cases.  The  pressure  penalty  increases  by 
about  10%  for  the  33%  GDL  intrusion.  Interestingly,  the  predicted 
pressure  drop  for  33%  intrusion  is  closest  to  the  experimental  mea¬ 
surement.  The  experimental  pressure  drop  for  this  case  is  460  Pa 
and  numerically  predicted  pressure  drop  is  423  Pa.  At  the  same  time 
the  pressure  drop  for  perfect  channels  is  much  less  than  the  exper¬ 
imental  data.  This  suggests  the  presence  of  GDL  intrusion  in  flow 
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Fig.  11.  Liquid  water  volume  fraction  at  the  section  of  anode  channels. 


channels  and  the  extent  of  area  maldistribution  likely  to  be  close  to 
33%.  The  liquid  volume  fraction  contours  in  the  cathode  channels 
for  the  same  cases  are  displayed  in  Fig.  8.  The  dry-to-wet  transition 
is  clearly  visible  in  all  cases.  With  GDL  intrusion,  the  maximum  liq¬ 
uid  volume  fraction  occurs  in  the  intruded  channels  and  the  liquid 
front  is  pushed  considerably  upstream.  This  is  because  the  intruded 
channels  feature  larger  flow  resistance  and  hence  lower  gas  veloc¬ 
ity,  thereby  accumulating  more  liquid  water.  The  presence  of  more 
liquid  water  further  increases  the  flow  resistance  and  reduces  the 
gas  flow  through  the  channel.  This  feedback  mechanism  results  in 
a  ‘U’  shape  of  flow  distribution  across  seven  channels. 

The  liquid  water  saturation  in  the  midsection  of  the  cathode  GDL 
for  these  three  cases  is  plotted  in  Fig.  9.  It  is  seen  that  the  liquid  sat¬ 
uration  is  much  higher  in  the  GDL,  reaching  about  30%.  Maximum 
liquid  water  saturation  appears  over  the  lands  at  about  the  mid¬ 
length  of  the  channels.  For  the  unintruded  channels  (Fig.  9a)  the 
liquid  water  saturation  profiles  exhibit  a  similar  pattern.  Interest¬ 
ingly,  Fig.  9b  and  c  show  that  the  liquid  water  saturation  in  the  GDL 
decreases  in  the  vicinity  of  the  intruded  channels,  possibly  due  to 
lower  current  density  (see  Fig.  6b  and  c).  The  effect  of  GDL  intrusion 
on  cell  operation  can  be  summarized  by  comparing  Figs.  7-9.  The 
current  production  is  low  over  the  intruded  channels  due  to  low 
reactant  flowrate.  As  a  result,  water  production  is  low  and  there¬ 
fore  the  liquid  saturation  in  the  GDL  is  low  over  the  intruded  edge 
channels.  Intuition  suggests  that  the  liquid  water  would  be  less  in 


the  intruded  channels  since  less  water  is  produced,  but  the  ability 
of  these  channels  to  remove  liquid  water  out  of  the  cell  is  so  low 
that  more  liquid  water  accumulates  in  the  intruded  channels. 

The  pressure  contours  along  the  anode  channels  are  shown  in 
Fig.  10  for  the  same  three  cases.  Pressure  penalty  also  increases  in 
the  anode  side  due  to  the  presence  of  two-phase  flow.  However, 
the  magnitude  of  the  anode  pressure  drop  is  about  20%  of  the  cath¬ 
ode  due  to  much  lower  hydrogen  flowrate  than  air.  The  liquid  water 
volume  fraction  contours  along  the  anode  channels  are  displayed  in 
Fig.  11.  Not  only  is  the  dry-to-wet  transition  captured  in  the  anode 
channels  but  there  is  also  a  wet-to-dry  transition.  The  latter  tran¬ 
sition  occurs  because  the  anode  outlet  is  facing  the  cathode  inlet 
where  relatively  dry  air  is  fed  into  the  cell  which  re-evaporates  the 
liquid  water  in  the  anode  channels.  It  is  seen  that  GDL  intrusion  has 
a  severe  effect  in  the  anode. 

The  liquid  water  saturation  in  the  midsection  of  the  anode  GDL 
for  these  three  cases  is  plotted  in  Fig.  12.  It  is  seen  that  the  liquid 
saturation  is  much  higher  in  the  GDL  than  the  anode  channels.  The 
maximum  liquid  saturation  occurs  over  the  lands  in  the  middle  of 
the  cell,  ranging  from  23%  to  25%.  Once  again,  for  the  unintruded 
channels  (Fig.  12a)  the  liquid  saturation  profiles  are  relatively  sym¬ 
metric  from  channel  to  channel,  while  there  is  clearly  more  liquid 
water  in  the  GDL  over  the  land  area  towards  the  center  of  the  cell 
with  intruded  edge  channels.  This  behavior  of  liquid  water  maldis¬ 
tribution  resembles  the  cathode  side. 
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Fig.  12.  Liquid  water  saturation  contours  at  the  midsection  of  the  anode  GDL. 


Flow  and  liquid  water  maldistribution  can  be  better  depicted  by 
liquid  water  volume  fraction  contours  in  the  center  of  the  channels 
along  the  in-plane  direction.  Fig.  13  compares  these  contour  plots 
for  a  middle  channel  and  an  edge  channel  for  the  three  cases  of 
interest.  For  the  perfect  channels  the  axial  profile  of  liquid  water  is 
very  similar  between  the  two  channels  (Fig.  13a).  Along  the  cathode 
channel,  a  dry-to-wet  transition  can  be  seen,  while  the  cathode  exit 
always  remains  wet.  In  the  anode  channels,  however,  there  exists 
dry-wet-dry  transition,  the  physics  of  which  has  been  elaborately 
discussed  in  the  literature  [32,35],  albeit  in  the  absence  of  two- 
phase  flow  in  channels.  Our  results  illustrate  a  profound  effect  of 
two-phase  flow  in  channels  on  the  liquid  water  distribution  in  the 
PEFC.  That  is,  the  wet-to-dry  transition  in  the  cathode  is  predicted 
without  channel  two-phase  flow  but  not  with  the  consideration 
of  channel  two-phase  flow.  This  is  because  the  channel  two-phase 
flow  allows  for  the  accumulation  of  a  large  amount  of  liquid  water, 
which  then  becomes  difficult  to  dry  up  by  dry  gas  from  the  anode 
inlet.  For  the  same  reason,  however,  the  anode  flow  stream  is  prone 
to  dryout  because  the  anode  flowrate  is  only  about  one-fifth  of  that 
of  the  cathode  side. 

The  effect  of  area  maldistribution  on  liquid  water  distribution 
is  more  evident  from  Fig.  13b.  It  is  shown  that  the  amount  of  liq¬ 
uid  water  is  much  higher  in  the  edge  channel  section  for  17%  GDL 
intrusion.  The  dry-to-wet  transition  in  the  cathode  occurs  further 
upstream  in  the  edge  channel  with  respect  to  the  middle  chan¬ 
nel.  For  the  case  with  33%  GDL  intrusion  these  effects  are  further 


magnified  as  shown  in  Fig.  13c.  The  difference  in  the  maximum 
liquid  water  volume  fraction  between  the  channels  is  quite  sub¬ 
stantial.  The  wetted  area  is  much  higher  for  edge  channels  as  well. 
Schematic  diagrams  illustrating  the  dramatic  effect  of  GDL  intru¬ 
sion  on  liquid  water  distribution  are  given  in  Fig.  14a  and  b  for 
the  middle  (unintruded)  and  edge  (intruded)  channels,  respec¬ 
tively. 

The  local  stoichiometry  in  each  of  the  seven  channels,  defined 
as  £avg(Q/Qavg)  with  Q. being  the  local  gas  flowrate  through  a  certain 
channel,  is  shown  in  Fig.  15  to  quantify  the  degree  of  flow  maldis¬ 
tribution.  Such  plots  for  both  cathode  and  anode  sides  are  shown 
in  Fig.  15.  The  channels  are  numbered  starting  from  the  inlet  in 
each  case.  It  can  be  seen  that  the  flow  maldistribution  is  relatively 
minor  in  perfect  channels  and  arises  solely  from  two-phase  flow. 
When  GDL  intrusion  is  considered,  severe  maldistribution  of  flow 
develops  in  the  edge  channels.  Area  reduction  of  17%  at  the  edge 
channels  results  in  about  25%  reduction  in  flowrate.  Similarly,  for 
33%  area  reduction  at  the  edge  channels,  the  reactant  flow  decreases 
by  almost  50%,  while  the  flows  through  the  other  channels  with¬ 
out  GDL  intrusion  do  not  vary  much.  These  observations  are  true  for 
both  the  cathode  and  anode  channels.  In  the  present  cases,  the  aver¬ 
age  flow  stoichiometry  over  the  entire  cell  is  2.0  in  both  anode  and 
cathode.  A  reduction  of  50%  flow  means  that  the  local  stoichiome¬ 
try  is  less  than  1.0,  insufficient  to  sustain  the  average  reaction  rate 
and  potentially  leading  to  detrimental  side  reactions  such  as  oxygen 
evolution  and  carbon  corrosion. 
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Fig.  13.  Liquid  water  saturation  contours  across  edge  and  middle  channels. 
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Fig.  14.  Schematic  of  liquid  water  distribution  in  a  PEFC. 
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Fig.  15.  Plots  of  local  flow  stoichiometry  in  different  channels  on  the  cathode  and 
anode  sites. 

6.  Conclusion 

A  channel  two-phase  flow  model  has  been  integrated  with  the 
previously  developed  two-phase  PEFC  model  based  on  the  M2 
framework.  The  complete  PEFC  model  considering  two-phase  flow 
in  channels  is  employed  to  explore  flow  maldistribution  in  an  oper¬ 
ating  PEFC  for  the  first  time.  The  wetted  area  ratio  on  the  cathode 
GDL  surface  predicted  by  the  present  model  matches  quantitatively 
with  experimental  data  over  a  range  of  current  density,  relative 
humidity  and  flow  stoichiometry.  In  addition,  the  overall  pressure 
drop  prediction  is  found  to  be  good. 


The  effect  of  GDL  intrusion  at  the  edge  channels,  as  commonly 
observed  in  a  fuel  cell  stack,  is  numerically  studied.  Severe  flow 
and  liquid  water  maldistributions  are  predicted  due  to  GDL  intru¬ 
sion  in  the  edge  channels.  Low  flowrate  of  the  intruded  channels 
make  these  regions  starved  of  reactants,  thus  reducing  the  cell  per¬ 
formance  and  durability.  Innovative  flow  field  designs  are  needed 
to  mitigate  flow  maldistribution  and  ensuing  adverse  impact  on  cell 
performance  and  durability. 
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